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Abstract 

Based on recently derived exact stochastic Liouville-von Neumann equations, several 
strategies for the efficient simulation of open quantum systems are developed and 
tested on the spin-boson model. The accuracy and efficiency of these simulations 
is verified for several test cases including both coherent and incoherent dynamics, 
involving timescales differing by several orders of magnitude. Using simulations 
with a time-dependent field, the time evolution of coherences in the reduced density 
matrix is investigated. Even in the case of weak damping, pronounced preparation 
effects are found. These indicate hidden coherence in the interacting system which 
can only be indirectly observed in the basis of the reduced quantum dynamics. 
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1 Introduction 



The model of a two-level atom interacting with a quantum field or many- 
particle reservoir has shown a remarkable persistence over decades of progress 
in physics, appearing and reappearing in many guises in different branches of 
condensed-matter physics after its 'first life' in quantum optics and magnetic 
resonance. It constitutes an idealized, minimal model of quantum dynamics 
and thermodynamics of open quantum systems. Through the explicit inclusion 
of a reservoir in the model, it supports the discussion of an open quantum 
system without resorting to speculative extensions of quantum mechanics. 

The spin-boson problem [1,2,3,4] is fairly well understood for most parameter 
regimes relevant to solid-state physics, where the reservoir is characterized by 
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a smooth spectral density running up to a large bandwidth cutoff, which en- 
ters the effective physics only as a renormalization parameter. A rich body of 
theoretical work on the spin-boson problem was developed in the context of 
macroscopic quantum coherence [5] and its realization in superconducting de- 
vices as well as the physics of defect tunneling in solids [6], to name two early 
examples. These results are currently being reapplied and extended in the dis- 
cussion of performance limits of quantum computers, since the two-level atom 
of the spin-boson problem can obviously be identified with a qubit subject to 
a dephasing-inducing environment [7]. Apart from this general consideration, 
the spin-boson model also serves as a model for particular realizations, e.g., 
the flux qubit [8], where physical parameters can be given for the two- level 
system and its environment. Lastly, it should be mentioned that an intricate 
formal link has recently been pointed out [9,10] between the dynamics of the 
spin-boson model and both dissipative transport of light particles [11] in a pe- 
riodic potential and transport of correlated electrons through a barrier in a 
ID conductor [12]. 

In the context of chemical physics, charge transfer and curve-crossing problems 
form an important context for the spin-boson problem. Here the understand- 
ing of the spin-boson dynamics seems less complete: The bandwidth of the 
dissipative reservoir may be small, or its spectrum may exhibit structured 
features relating to vibrational spectra. In addition to the energy scales of the 
two- level system and a damping constant, other parameters need to be con- 
sidered. Moreover, the spin-boson treatment of chemical phenomena is often 
considered only a first approximation, since the linear-response assumption 
for the reservoir is inherent in the model; a generalization to nonlinear, more 
complex reservoirs is sometimes warranted. Since there is no mature theory of 
non-linear reservoirs, the straightforward way to achieve such a generalization 
is the formal inclusion of key degrees of freedom of the reservoir, e.g., some 
prominent vibrational modes, in the open system. But even the spin-boson 
problem itself is not fully understood absent the scaling behaviour found in 
solid-state physics. 

It is not to be expected that a thorough theoretical treatment of the more 
complex models just indicated can be accomplished by analytic methods alone; 
numerical simulations will become more important as larger open systems are 
studied. Likewise, simulations are likely to help in charting the remaining terra 
incognita of the spin-boson problem, e.g., the electron transfer problem in the 
inverted regime. Future progress will place competing demands on numerical 
methods which most established algorithms cannot meet at the same time: 
acceptable scaling of computational complexity with both system size and time 
interval simulated [13]. Additionally, the method should take into account the 
generally non-Markovian nature of quantum fluctuations. 

The usual starting point for the theoretical description of an open quantum 
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system beyond perturbation theory, by now more or less canonical, uses the 
path integral formalism, which allows the exact treatment of quantum memory 
effects [14]. Since quantum memory effects do not allow the path integral to 
be translated into a simple differential equation of motion, this limits the 
choice of theoretical and numerical approaches available. The direct evaluation 
of path integrals using Monte Carlo sampling is a reliable method for short 
to intermediate times [15,16,17] but becomes prohibitively expensive at long 
times due to the dynamical sign problem. The recursive evaluation of a path 
integral using the quasiadiabatic path-integral discretization (QUAPI) allows 
propagation to arbitrarily long times, but needs computational resources which 
grow very rapidly with increasing system size or memory time [18,19,20]. 

More recently, stochastic approaches have allowed a transition from the path 
integral description to equations of motion [21,22,23,24,25,26]. Similar to the 
intuitive - and computationally advantageous - representation of classical dis- 
sipative systems through the motion of a phase-space point under the influence 
of thermal noise and friction, an open quantum system can be described by 
the stochastic propagation of pure quantum states [21,22,23,24,25,26,27,28]. 
Favourable scaling with system size seems an intrinsic feature common to 
these approaches - the computational complexity of the open-system simula- 
tion as a function of system size remains the same as that of simulating the 
same quantum system without external interaction. The remainder of this ar- 
ticle is largely devoted to the question how the second objective, good scaling 
for long-time dynamics, can be achieved in the presence of memory effects. For 
tests of this performance aspect, the spin-boson system is a prime candidate, 
since comparison results are readily available. 

A summary of the formal description of dissipative quantum systems through 
stochastic Liouville-von Neumann (SLN) equations given in Refs. [25,26] in 
Section 2. Section 3 discusses new material, giving details of numerical ap- 
proaches based on the SLN formalism. A series of numerical tests on the free 
dynamics of the spin-boson model as well as simulation results probing hidden 
coherence in the spin-boson model through pulsed external fields are presented 
in Section 4, followed by a summary of results and conclusions. An appendix 
presents the derivation of a mathematical result needed in Section 2. 



2 Stochastic Liouville-von Neumann equations for open quantum 
systems 

Although a connection between the influence functional formalism and classi- 
cal coloured noise in quantum dynamics was pointed out already in the seminal 
work of Feynman and Vernon [14], it has been put to use in the theory and 
in simulations of open quantum systems only recently [21,22,23,24,25,26] (see, 
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however, [29,30]). The key ideas of one such stochastic approach to open quan- 
tum systems [25,26] will be outlined in this section using a more general model 
than the spin-boson model, commonly known as the Caldeira-Legget model. 
The Caldeira-Legget model consists of a one-dimensional potential model cou- 
pled to a thermal reservoir of harmonic oscillators with a quasicontinuous 
distribution of frequencies, 

2m w Y 2 V m jUj) 2m i 



With the identification 

1 ~> 7T^' ( 2 ) 



where q is a characteristic length scale, the spin-boson Hamiltonian 

hA he go V- rrij-Wj 2 

#sb = - + —a z - —<r z ^ cj-xj + ^ + ~2 ^ ( 3 ) 

3 3 3 



can be considered a truncation of the Caldeira-Legget model to two localized 
states [31,3]. 

The dynamics of the harmonic reservoir is usually of little interest; moreover, 
due to the large size of the total system, only a reduced description from 
which the oscillators have been eliminated is simple enough to be practically 
treatable. Although Master equations have been frequently been used with 
considerable success in weak-coupling scenarios such as quantum optics [32], 
a formally exact reduced formalism, usually necessary in a condensed-matter 
context, is known only in path integral formalism. Feynman and Vernon [14] 
have demonstrated that the path integral 



r,' 

q t If 

p(qt,&t) = jdq, jdq[jD[ qi ] jD[q 2 ] 

9i q[ 

xexp Q(Sb[gi] - Sb[g 2 ])) 

xF[ qi - q 2 , ( Ql + q 2 )/2] p(q u q[; t ) (4) 

is an exact expression for the time evaluation of the reduced density matrix 
p{q{,q'f',t) for an initial density matrix which factorizes between p(<fc, to) 
and a thermal density matrix of the harmonic reservoir. Here So[q] denotes 
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the classical action functional associated with the one-dimensional potential 
model. The influence functional 



F[x, r] = exp I — — r] 



(5) 



with 




dt' dt"x(t') 



x [Re L{t' - t")x(t") + 2i Im L(t' - t")r(t")} 



t 



+-£ dt' x(t')r(t') . 




(6) 



contains the full physics of the system-reservoir interaction and all aspects of 
the reservoir dynamics which may be reflected in the system dynamics. The 
kernel 



describes free fluctuations of the reservoir's coupling coordinate, themselves 
dependent on a spectral density J{oj) and inverse thermal energy j3. Typical 
spectral densities in the context of solid-state physics are smooth in the rel- 
evant frequency regime, e.g., proportional to u 3 for bulk phonons or of the 
Ohmic form J(oo) = rjuj for the low-energy excitations of a Fermi liquid. In 
the case of the spin-boson model, a dimensionless constant a = q^rj/^nh) 
describes the strength of Ohmic dissipation. 

The last term in Eq. (6) has the form of a potential modification, which re- 
flects the fact that the system-reservoir interaction given in Eq.(l) eliminates 
any quasistatic response of the reservoir; it relates to dynamic reservoir's 
dynamic response function xn{t ~ t') = —2Q(t — t') \m.L(t — t')/% through 
\i = J °° drxR(r). This potential term has no effect for the spin-boson model, 
from q 2 = q^/A one finds fixr = (n/2)(q 2 — q%) = 0. 

The path integral expression (4) has the distinct advantage of providing an 
exact description of the system-plus-reservoir dynamics which is reduced to 
only the system degree of freedom. However, the propagator for the reduced 
density matrix described by Eq. (4) is not associative due to the fact that 
the exponent of F[x,r] contains a double time integral, i.e., there are memory 
effects which make quantum amplitudes a non-local functional on the path 
space. 



L(t) = - [du J(cj)(coth — 

7T J 2 



cos out — i sin cut) 



(7) 
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At the price of introducing a further functional integral over auxiliary function 
spaces of complex functions £(t) and u(t), this non- locality can be lifted, repre- 
senting the influence functional as a weighted average of quantum amplitudes 
containing a time-local action functional 



F[x,r] = J D 2 l$ J I?[u]W[Z,?,u,u*] (8) 
xexpl^J dt'Z(t')x(t') + \u{t')r{t' 

x exp J dt'x(t')r(t') 



to 



with a suitable Gaussian functional £*, v, v*\. Because W is normalized, it 
may be interpreted as a probability density of real noise fluctuations (£ +£*)/2, 
(£ _ C*)/(2i), + ^*)/2 and (1/ - v*)/(2\). It is partly characterized by the 
noise correlation functions 



(atW)) w = ReL(t-f) (9) 
(Z(t)v(t')) w = (2i/h)e(t - t') lmL(t - t') 

= -i XR (t-t') (10) 

(u(t)u(t')) w = (11) 

required for the identification of Eqs. (8) and (5). W is also characterized by 
the correlations (£(t)£*(t')), (£(t)v*(t')) and (v(t)v*(f)), which do not enter 
the physical result obtained upon stochastic averaging. 

The immediate benefit of the stochastic construction (8) becomes clear if the 
order of the integrations over (q, q') and (£, v) is interchanged: for any specific 
functions £(t) and u(t) the path-integral dynamics can be translated into the 
Schrodinger picture in the usual way. Observing that the exponents in Eq. 
(8) can be identified as action functionals associated with a time-dependent 
potential, we immediately recover the equation of motion 

ihp = [H 0: p}_ - £(*)[?, d- + \\q\p\- - ^(t)[q,p} + , (12) 



a stochastic Liouville-von Neumann (SLN) equation. The presence of both 
complex parameters and an anticommutator makes this equation describe a 
non-unitary propagation of individual samples, allowing p to acquire a non- 
hermitean component. After stochastic averaging, however, any non-hermitean 
parts of p vanish. Similarly, the trace of p varies between individual samples, 
but remains equal to unity on average. 
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This diffusive spread of the sample trace can be reduced or eliminated by 
introducing a 'guide trajectory' (to be specified later) in the anticommutator 
term, 

ih'p = [H ,p]--S(t)[q,p]- + ^{q\f>V - ^u(t)[q-f t ,p]+, (13) 



where 



p(t) = explij dt'v(t')f t ' p(t) (14) 



relates p to the solution of the original equation (12) for arbitrary f t . It is 
to be noted that Eqs. (13) and (14) do not set the stage for an approximate 
expansion; together they form an exact identity. The guide trajectory may 
depend on the noise forces in the form of a functional f t [C,(t'), v(t')], where the 
following properties will be assumed: 



() \ = 0;if>t (15) 



8£(t 
Sr t 



M*0 = ° ;t '-' (16) 

= 0; any t'. (18) 

In words, the functional r t [£(i'), v(t')] is required to be causal and analytic. 

Under these conditions, the exponential factor in Eq. (14) can be absorbed 
into the probability measure by re-defining the centre of the Gaussian func- 
tional W[£, £*, v, v*\ (see Appendix). This leads to new noise forces with same 
variance but dynamically shifted mean values, 



o 

t 

= £(0 + jds XR (t' -s)f(s) (19) 

o 
t 

C t (t') = C(t')-^Jds(C(t'HsMs) (20) 
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vt{t') = v{t') (21) 

t 

v * t (t') = v*{t') - i jds (v*{t')v{s))r{s) (22) 
o 

Dropping the subscript t after the transformation has been performed, we 
obtain the SLN equation 



ih'p = [H , p]_ - (t(t) - J dt' XR (t - t')f t )j [q, p]_ (23) 

to be averaged with the original probability density W[£, £*, v, v*\ without the 
exponential factor of Eq. (14). 

A natural choice for f t relates it to the dynamics on the time interval [0,t], 



Since this makes Eq. (23) conserve the trace of p for each sample, we call 
Eq. (23) a normalized SLN equation for a guide trajectory defined by Eq. 
(24). Compared to the linear SLN equation (12) the normalized version allows 
a more efficient stochastic averaging as well as a clear interpretation of the 
classical limit of Eq. (23) [25,26]. 

The implementation of a numerical simulation algorithm based on the normal- 
ized SLN equation is now relatively straightforward. For a numerical construc- 
tion of the noise fluctuations, £ is decomposed into a purely real £^ with a 
long autocorrelation timescale, which is uncorrelated to all other noise forces, 
and a short-range complex part t^ s \ with correlations 



(^\t)^(t')) = ReL(t-t') (25) 

(e W (f)£ W (O> = (26) 

(Z (s) (tnti)) = -i X R(t-t>) (27) 

(£«(fK(f)> = (28) 

(v(t)v(t)) = 0. (29) 



Here the correlations of ^ and v decay rapidly for uj c \t — t'\ ^> 1, while 
correlations of £®(t) persist up to the thermal timescale h{3. 
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These noise fluctuations are generated numerically with great efficiency by fil- 
tering white noise with suitable integral operators, whose kernels have spectra 
which multiply to yield the spectra of L(t — t') and Xn{t ~ t')- The neces- 
sary convolution operations can be performed with ease using a fast Fourier 
transform algorithm. This construction also determines the non-physical cor- 
relations (£ (s) (t)£ (s) *(t')> and (i/(f)i/*(f)). 

The nonlinearity of the numerically simulated dynamics can be reduced by 
observing that the solution of the simplified SLN equation 



\Kp=[H ,p\- - [i{t) -J dt' XR (t - t')f t )j [q,p] 



(30) 



immediately yields a solution of the nonlinear SLN equation (23) through the 
relation 

m - s^/k*). (3D 



A further simplification of Eq. (30) lies in the factorizing ansatz p = IV'iKV^Ij 
which reduces Eq. (30) to the two Schrodinger equations 



t 



ih\ij> 1 ) = Ho\1> 1 ) - \ t(t) -J dt' X R(t-t')r t j q\^) 



+^\^} - \v{t)q\^) (32) 



t 



ih\i, 2 ) = H,\i, 2 ) - \ ?{t) -J dt' XR (t-t')f* t ,j q\ij 2 ) 



+ f ?V2> + \s{t)q\fr) (33) 



Apart from the satisfaction of arriving at a conceptually simple result from an 
involved mathematical transformation, this reflects very positively on the scal- 
ing of stochastic simulations: For a single noise sample, the numerical effort 
of propagating the system scales with system size exactly the same way as for 
a conservative system. Moreover, it is interesting to note that Eqs. (32) and 
(33) capture all effects of the system-reservoir interaction, including quantum 
correlations and entanglement between system and reservoir, through the ad- 
dition of single-particle operators to the Hamiltonian H . The description of 
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quantum dissipation through single-particle operators strongly suggests that 
a wide range of approximations to conservative quantum systems (e.g., for- 
malisms suitable for large atoms or small molecules) can be combined with this 
stochastic approach. Larger and more complex quantum systems, which can 
currently be simulated only under the assumption of energy conservation, are 
likely to become accessible to an open-system approach using SLN equations. 
There is an attractive feature in in the SLN approach which may recommend 
it for complex systems even in the case of very weak damping, where perturba- 
tive methods are valid: Whereas Born-Markov perturbative approaches such 
as Redfield equations define damping in terms of system and reservoir char- 
acteristics, often making reference to an exact diagonalization of the system, 
the SLN noise terms are constructed from the free reservoir dynamics alone; 
no characteristics of the free system dynamics enter the formal description of 
dissipation. 



3 Optimized stochastic simulation methods 



For the transformed noise forces given in Eqs.(15) - (18) the linear combi- 
nations £t ± £j and v t ± are no longer purely real or purely imaginary. 
The transformed integration measure is formally Gaussian, but the integra- 
tion, even when written over the real components of £ and z/, is shifted in the 
complex domain. However, the stochastic simulation of Eq. (30) must use real- 
valued components for these shifted quantities, i.e., the numerical simulation 
of Eq. (30) operates on an analytic continuation of the exact result. Although 
this continuation is known to be free of singularities in important cases such 
as the classical limit, the weak-coupling limit and any harmonic system, the 
spin-boson model (at low temperature) appears to be a counterexample where 
the simulation becomes unstable after returning excellent results up to a well- 
defined time threshold t c , which is related to the growth of the non-hermitean 
part of p. Beyond the threshold a sudden onset of peak-like artefacts and other 
systematic errors is observed, in close similarity to stability problems found 
in the Master equations for the positive P function, a quantum optical phase 
space function [33,34]. As in the latter example, SLN simulations using Eq. 
(30) are highly accurate for times below the threshold time. In the light of 
these facts two natural ways to extend the power of the SLN methodology are 
(a) increasing the timescale t c and (b) finding consistent methods to 'reshuffle' 
the ensemble to a healthy state before artefacts taint the data. 
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3.1 Alternative noise spectra 



In the scaling limit any relevant frequency Vt of the dissipative dynamics lies 
within the band of reservoir frequencies, Q <C u c , i.e., the complex noise 
forces, which drive the system away from normalized hermitean states, act 
on the system at resonance. This can be changed if the potential term in the 
influence functional is formally included in the stochastic construction, which 
leads to modified cross-correlations for the noise fluctuations ^ and v, 

(&\t)v{t>)) = - 1XR (t - t') + itf(t - 0- (34) 



Instead of a noise spectrum extending roughly over the interval [— u c ,ui c ], we 
are now dealing with a white noise spectrum that has a gap in the same 
interval. The resonant driving of the system by complex noise is thus strongly 
reduced, typically increasing the critical time t c by one or more orders of 
magnitude. With the proviso that increments of the white noise parts of £ and 
v be interpreted in the sense of a Stratonovich stochastic integral, the steps 
taken in the derivation of Eqs. (23) and (30) can be reiterated unchanged. 
This yields 



ihp = [H , p}_- \ £(t) - J dt' XR (t - t')f t , + fj,r t j [q, p]_ 

~u(t)[q,p} + . (35) 

as a modified version of the simplified SLN equation (30), which has the ad- 
ditional advantage of having only translationally invariant dissipation terms. 



3.2 Collective normalization of subensembles 



Next we consider the case of a subensemble of iV samples p^ (t), each governed 
by an independent set of noise forces v^\t). In this case, a set of 

guide trajectories can be given which does not conserve the trace of each 
sample, but only their sum. Again, guide trajectories enter an exponential 
factor distinguishing normalized and unnormalized states, but here a single 
factor is used for different noise realizations, 

N 

p U)(t) = J] exp ( i / dt'^ k \t')f ( t k) ] p( j \t') (36) 
k=i 
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with 



M _ trgpW(t) _ 1 ^ (fc) 
* "E^trp(0(t) TV 



trgpW(t) (37) 



For initially normalized states, this keeps the sum over traces identical to N. 

Again, the exponential factors in Eq. (36) must be absorbed into the proba- 
bility measure through suitable shifts of the noise fluctuations. Because the 
(unshifted) noises acting on different members of the subensemble are uncor- 
rected, and because terms with different in Eq. (36) factorize, the shifts 
of noise fluctuations v®, v®*) and (£(*),£(*)*, i/ fc ), i/*)*) with j ^ k 

are independent of each other. Thus the derivation presented in the Appendix 
applies also in this case. Eq. (19) is only trivially modified for collectively 
normalized samples, 

t 

tf\t>) = ) + Jds XR {t' - s)f?{s). (38) 



The corresponding normalized SLN equation reads 



tip® = [H ,p®]- - ^ 0) (0 - / dt'xn{t - t')ff^ [q,p®\_ 

+^k 2 ,p]- - ^ U) (t)hp U) ] + + E h^KttfVpV. (39) 

k=l 

Again, there is a simplified SLN equation suitable for numerical simulation, 



ift/iW = [H ,p^}_ - U<fi(t) - J dt' XR (t - t'rfp] [q,p^U 



to 



+f[? 2 ,P]--^ (, ' ) (*)[?,P (, ' ) ] + , (40) 
related to the normalized form (39) by 

= ^Jp^/ ' {t) - (41) 



It is evident from Eq. (37) that is of order 1/N, i.e., the shift (and the 
corresponding nonlinearity of the equation of motion) vanishes for large N; 
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Fig. 1. Schematic representation of block factorization of short-time correlations 
(left) vs. global reach of long-time correlations (right) (see subsection 3.3). 

the threshold time t c is thus further extended. The linear SLN equation (12) 
is recovered from (40) in the limit iV — > oo. 

3.3 Separation of friction and thermal time scales 

The split of £ into the long- and short-range components £W and £w allows the 
numerical solution of SLN equations using a two-stage averaging procedure, 
which separates the difficulty of dealing with long-range correlations from the 
problems incurred through complex driving forces. In a primary averaging 
stage (inner sampling loop), samples are drawn for £W and u, while £® will 
be changed only in the secondary averaging stage (outer sampling loop). For 
simulation parameters yielding a sufficiently large threshold time, u c t c 1, 
the short-time noise correlation (£,(t)v(t')) can be block factorized on intervals 
of width r < t c , as indicated in Fig. 1 (left). Under the condition uj c t ^> 1 this 
leaves the effective support of the correlator (£(t)z/(t')) virtually untouched. 

In the primary averaging, the dynamical simulation can thus be partitioned 
over finite time intervals, accumulating and then drawing new Hermitean sam- 
ples at the end of each interval. This procedure does not disturb the long-time 
correlations since all samples drawn for the primary averaging evolve subject 
to the same realization of In the secondary averaging, the entire procedure 
is repeated for a sufficient number of long-range fluctuations which remain 
unconstrained in the t,t' plane. Their effective support forms a diagonal band 
of width oc h/3 (Fig. 1, right). 



4 Tests and results for the spin-boson model 

Numerical tests and applications of the simulation strategies described in Sec- 
tion 3 are discussed in the following; a combination of all three strategies is 
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used unless indicated otherwise. The underlying SLN equations of motion for 
the spin-boson model are as follows: The linear SLN equation reads 



, . hA he, , J _. . r . hv 

9 = — ^ ax,p >~ + y^'d- ~ twl?*' p\- ~ yp*-^ ( 42 ) 



Absent any physical context defining a characteristic length g , is identified 
with the position operator q. The normalized and simplified SLN equations 
are 



+ ~ hA . ^ he r ^ 



to 

— 7rWz-<Tt,p\+ (43) 



2 

(with <7 t = tr o z p) and 



ftA . . foe r . 



- [i{t) -J dt' XR {t-t')a)j [a z ,p}_ 



(44) 

(with p — p/trp), respectively. Throughout the following, Ohmic dissipation 
with an algebraic cutoff 



+ 



corresponding to an exponentially decaying response function 

XR(t) = /j,Uctexp(-uj c t) (46) 



is used. 
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4-1 Coherent dynamics: free oscillations 

The coherent dynamics of the spin-boson model at moderate damping, a < 
1/2, and low temperature are theoretically well understood and numerically 
accessible using other methods, thus providing a suitable test case to verify 
the accuracy and stability of the SLN approach. Using data obtained through 
an established path- integral Monte Carlo (PIMC) method for comparison [15], 
the power and efficiency of the SLN approach in the region of coherent dynam- 
ics at zero temperature is demonstrated. Fig. 2 shows the time evolution of 
the expectation value (a z ) of a symmetric system from an initially factorized 
state. 

In the SLN approach the dynamics can be accurately simulated over many 
oscillation periods; empirical data suggest that the growth of statistical er- 
rors saturates to a constant value for increasing t, allowing the simulation to 
be continued to arbitrarily long times. This overcomes an inherent restric- 
tion of the PIMC approach, which suffers from the so-called dynamical sign 
problem: The PIMC approach mainly employs the superposition principle to 
obtain quantum amplitudes. The resulting averaging over highly oscillatory 
functions results in a signal-to-noise ration of the MC simulation which de- 
grades exponentially with growing t, leading to near-divergent error estimates 
after about one oscillation period. For very strong system-reservoir coupling, 
however, quantum phase factors are suppressed, and the PIMC approach gains 
in efficiency. 

The computer times of the SLN and PIMC simulations shown in Fig. 2 are 
comparable; they vary between 40 h and 210 h on a commodity CPU (Athlon, 
1200 MHz). The SLN simulations, whose statistical errors are comparable to 
the linewidth, use between 10 and 20 factorization blocks. The time axis is 
scaled by the renormalized matrix element A r = A(A/u c ) ai/( - 1 ~ a \ In the SLN 
simulations a high cutoff uo c = 20 A was used while a lower cutoff u c = 6 A was 
used in the PIMC simulation to avoid ergodicity problems. 



4-2 Incoherent relaxation: long-time dynamics 

Simulating the slow exponential decay of localized populations in a biased 
spin-boson model, as it occurs, e.g., in electron transfer reactions, is chosen as 
another test case. For parameter regions with non-trivial transient behaviour, 
an explicit simulation of time-dependent population numbers can offer advan- 
tages over simulations linking rates to thermal flux correlation functions [17]. 
The population difference P(t) is chosen as a dynamical observable; P(t) is 
formally defined as the time-dependent expectation value o~ z (t) for a two-state 
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Fig. 2. Comparison of SLN simulation data (lines) and path-integral Monte Carlo 
data (symbols) for coherent spin-boson dynamics at zero temperature. The time t 
is scaled by the renormalized matrix element A r (see text). 

system starting from the factorized initial state with 



and the reservoir in thermal equilibrium. For rate calculations, it is essential 
that P(t) can be simulated over a long enough time interval to be directly 
fitted to an exponential function. Since Pit) decays towards a nonzero value 
in a biased system, it is advantageous to symmetrize it under a reversal of the 
bias sign, e — > —e, i.e., study 



Because changing the sign of e is the same as interchanging the eigenstates of 
a z , this can be expressed as a single expectation value with an antisymmetric 
initial state 



Using this equivalent definition of the symmetrized expectation value P s {t) as 
the basis of numerical simulations, the data presented in Fig. 3 is obtained. 




(47) 



Ps(t) = -(P £ (t)+P-e(t)). 





(49) 
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Fig. 3. Exponential relaxation of the asymmetric spin-boson system. All curves agree 
with analytic nonadiabatic rates within a few percent. The time t is given in units 
of A. 

Parameters are a = 0.1, u c /A = 100, k^T/TiA = 10 and e as indicated in the 
figure. Since the simulation time is many orders of magnitude longer than the 
inverse of u c , the relatively large number of 300 factorization blocks may be 
used. In spite of the large number of time steps, the computation of a single 
curve in Fig. 3 took only about 15 min. 

For rates varying over more than an order of magnitude, a remarkably sta- 
ble simulation result for exponential decay is found, extending over many 
time constants. Given the fact that these curves result from the averaging of 
coherently propagated samples, this is truly remarkable. In the light of this 
successful example of a highly efficient simulation, it is to be expected that 
more elaborate simulations using this technique will be powerful enough to 
generate benchmark data for electron transfer rates which can be used to 
check theoretical approximations. 

4-3 Rapid dephasing: short-time dynamics 

The dephasing of a coherent superposition of quantum states is one of the 
fastest processes in spin-boson dynamics. It is mainly governed by the in- 
teraction part of the Hamiltonian (3) [35]. After preparing the system in an 
eigenstate of a x , again with a factorizing initial condition, dephasing expresses 
itself through the decay of the off-diagonal elements of the density matrix, here 
probed by measuring (o~ x (t)). The results given in Fig. 4 are for different val- 
ues a = 5, 2, 1, 0.5, 0.2, 0.1 and 0.05 of the dissipation constant (different 
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Fig. 4. Rapid decoherence of a symmetric superposition of localized states with a 
factorized initial density matrix and uj c / A = 100, k-^T/UA = 20 and with a varying 
from 5 (symbol curve at left) to 0.05 (symbol '+', at top). The time t is given 
in units of A. 

symbols, from bottom left to top right). 

As outlined in Ref. [35], a transition from a rapid Gaussian decay towards a 
slower, more complicated decay is observed when varying the damping con- 
stant from strong to weak damping. Excellent agreement between the analytic 
strong-coupling result (solid line) and simulation (filled dots) is observed. The 
statistical errors are less than ±0.02. 

4-4 Preparation effects: probing hidden coherence 

In the case of moderate damping a < 1/2 and low enough temperature, the 
initial rapid decay of a coherent superposition leads to a relaxed state with 
some residual coherence, which can be probed by time-dependent external 
fields. As is commonly done in elementary treatments of magnetic resonance, 
the reduced density matrix of the two-level system can be parameterized by a 

— # 

polarization or spin vector M, 



position of the eigenstates of a z corresponds to a spin vector of unit length 
pointing in the x direction. Allowing the system to evolve from a factorized 




(50) 




Preparation of the two-state systems as a symmetric super- 
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Fig. 5. Dynamical simulation with a symmetric superposition as an initial state 
(factorizing with the reservoir state). The relaxed state after a rapid initial dephas- 
ing (inset) is probed by a pulsed field in z direction at t = 0, which rotates the 
spin vector by an angle of n/2. Subsequently the residual coherence of the relaxed 
state reveals itself through oscillatory dynamics. The rise of {a z ) above the relaxed 
value of {a x ) is a counterintuitive observation, indicating that the quantity {a x ) 
underestimates the coherence present in the interacting system. Solid lines indicate 
dynamics for a delta pulse at t = 0, dashed lines a finite- width pulse with equal area 
and shape as indicated in the figure. Parameters are A = 1, oj c = 250, a = 0.05 and 
k B T = 0. 



initial state (as in the previous subsection), the direction of M is conserved 
due to symmetry, while the reduction of its length indicates dephasing. In the 
case of weak damping, however, the initial dephasing saturates at a finite value 
of (o~ x ). As shown in Fig. 5, this residual coherence can be probed through a 
pulsed external field in the z direction, which changes the relative phase of the 
superposition, or, in other words, rotates the spin vector by an angle of ir/2, 
leaving it antiparallel to the y axis (without changing its norm). The follow- 
ing free dynamics displays damped Rabi oscillations resembling a precession 
of the spin vector. However, the norm \M\ = \J ((J x ) 2 + (cry) 2 + (cr z ) 2 of the 
spin vector does not fall monotonously, as one might expect for a dissipative 
system. After a quarter oscillation period M points in the z direction, with 
M z exceeding the relaxed value of \M\ before the onset of oscillations. This 
indicates hidden coherence in the interacting system, which does not become 
visible in the reduced density matrix. Results for a delta-shaped pulse (solid 
line) and a finite-width pulse (dashed line), whose shape is indicated at the 
bottom of Fig. 5, are virtually identical. 
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Fig. 6. Comparison of coherent oscillations starting from a relaxed delocalized state 
(solid lines, same as Fig. 5) and a factorized initial state at t = with identical 
reduced density matrix p (dashed lines). A rapid initial decay of (a y ) is observed for 
the factorized initial state, but not for the relaxed state, resulting in a significant 
difference in the amplitudes of coherent oscillations, while the expectation value 
{a x ) seems to be unaffected by the choice of initial preparation. Parameters are 
A = 1, u c = 250, a = 0.05 and k B T = 0. 

To test this interpretation in terms of hidden coherence, the free decay of 
(cr 2 ) and (a y ) is compared in Fig. 6 for the two different preparations of the 
initial state at t — 0. The rotated relaxed state (delta-pulse case of Fig. 5) 
is compared to the dynamics starting with a factorized initial state, where 
the reduced density matrix at t = is identical in both cases. A significant 
preparation effect is found: In the factorized case, coherence appears to be 
drastically reduced compared to the initial preparation derived from a relaxed 
state, supporting the notion of hidden coherence. While a rapid initial dephas- 
ing is again observed for the factorizing initial state, it appears to be absent 
for the relaxed preparation (see inset of Fig. 6). 

Signatures of hidden coherence can be observed without taking a symmetric 
state for a starting point. Fig. 7 shows a simulation starting from a con- 
ventional factorized preparation involving a localized state, which has been 
commonly used as a standard in the spin-boson literature. After a quarter- 
period of free oscillatory dynamics, a pulsed field in z direction rotates the 
spin vector to align it with the x direction, effectively stopping the precession 
of M. A second pulse later aligns the spin vector with the y axis, allowing 
the system to resume its oscillatory dynamics. Again, the maximum of (o~ z ) 
exceeds the previous norm of the spin vector in the delocalized state between 
the two pulses. 
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2 4 6 8 t 10 

Fig. 7. Two-pulse dynamics revealing hidden coherence: After the free precession 
of the spin vector starting from a conventional localized/factorizing initial state, a 
pulsed field applied at t = 3.2 rotates the spin vector from the x — y plane onto the x 
axis. At t = 4, a second pulse rotates the spin vector onto the y axis, allowing the free 
precession to resume. The amplitude of the following oscillations demonstrates that 
the value of (a x ) observed between pulses significantly underestimates coherence. 
Parameters are A = 1, a = 0.1, uj c /A = 1000 and k-gT = 0. 

In this last example, the separation of timescales between system and reservoir 
is so pronounced that the linear SLN equation (42) with the alternative noise 
spectrum (34) can be used for effective simulations. The simulations discussed 
in this subsection have also demonstrated that the SLN approach is capable of 
performing simulations whose characteristic timescales differ by several orders 
of magnitude. 



5 Conclusions and outlook 

A bundle of new numerical simulation strategies is added to the recently in- 
troduced stochastic Liouville-von Neumann equations for open quantum sys- 
tems. First, judicious use has been made of the freedom these equations allow 
in the choice of complex noise spectra. The resulting noise spectra largely 
avoid driving the system resonantly, leading to markedly improved numeri- 
cal performance. Additionally, the normalization of subensembles instead of 
individual samples improves the stability of the method in a crucial way. Fi- 
nally, use is made of the separation of timescales between dynamical response 
(inverse bandwidth) and the thermal timescale, yielding stable and efficient 
simulations for arbitrarily long times. 
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Using the spin-boson model as a test bed, these simulation methods have 
been validated as an efficient approach to the dynamics of open quantum 
systems. Several examples demonstrate its long-time stability as well as its 
suitability for problems with disparate timescales extending over several orders 
of magnitude. 

A series of simulation results using pulsed fields to extract information on 
coherent superpositions suggests that the off-diagonal elements of the reduced 
density matrix are not a suitable measure of quantum coherence in the spin- 
boson problem; there appears to be hidden coherence related to the system- 
reservoir interaction. 

Since the stochastic Liouville-von Neumann equations discussed here can be 
factorized into Schrodinger equations, they can be applied to larger systems 
without any complications other than those present already for conservative 
systems. The simple expressions found for the stochastic force operators sup- 
ports this prospect. 
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Appendix: Dynamical transformation of the probability measure 

In order to transform the linear dynamics of Eq. (12) into the nonlinear norm- 
conserving form (23) we need to find a set of functions £t(i')> v t{t'), 
tf(t') such that 

W[£ u C, v t , v\\ = W[£, T, v, v*\ exp [i J dfV(f)rV j (51) 

and show that the Jacobian determinant of the variable change to the new 
functions as independent variables of the noise averaging is unity. 
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Using vector notation z = (£, £*, */, i/*), the probability density W can be 
formally written as 

W%f>X] = ^ ex P v 4/^7 ^W-*V(0,) (52) 



where iV is a normalization constant. Here the superscript t denotes a trans- 
posed matrix (not a hermitean adjoint) because the inner product of the under- 
lying space of real functions must be used. The four-by-four matrix M(t' — t") 
is related to the noise correlation matrix (z t (t)z(t // )) through 



J drM(t' - T)(z\T)z(t")) = 5(t' - t"), 



(53) 



i.e., M(t' — t") is the inverse of the noise correlation matrix. With Eqs. (52) 
and (53) it is easy to show that Eqs. (19) - (22) 'complete the square' in the 
exponent of W. 

The additional task of computing the corresponding Jacobian only requires a 
fairly short calculation. Due to the properties required of the guide trajectory, 
one finds the simplification 
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(54) 



Here the second and fourth column simplify due to conditions (17) and (18) 
while the third row reflects the fact that v is not shifted. 

Now the definition (19) of £ t immediately leads to 



^ = S(t' - t") + J ds XR (t> - s] 



5r(s) 



(55) 



Because Xn(t' — s) is a causal response function and j0tj vanishes for t" > s 
by construction, the integrand is nonzero only for t" < i! . The integral term 
in Eq. (55) therefore does not enter the determinant, and we find 



J 



(56) 
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